Nonlinear current response of an isolated system of interacting fermions 
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Nonlinear real-time response of interacting particles is studied on the example of a one- 
dimensional tight-binding model of spinless fermions driven by electric field. Using equations of 
motion and numerical methods we show that for a non-integrable case at finite temperatures the 
major effect of nonlinearity can be taken into account within the linear response formalism extended 
by a renormalization of the kinetic energy due to the Joule heating. On the other hand, integrable 
systems show on constant driving a different universality with a damped oscillating current whereby 
the frequency is related but not equal to the Bloch oscillations. 
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Despite its importance for various branches of physics, 
the real-time response of quantum systems remains in 
many aspects an unexplored field. Only recently have 
the time-resolved measurements provided important in- 
formation on the nonequilibrium short-time dynamics 
of correlated bulk materials [l| nanostructures [2| and 
optical lattice systems [3[. In contrast to the develop- 
ments in the experimental methods, theoretical descrip- 
tion of the real-time dynamics remains a difficult and 
challenging task. As the exact time evolution is known 
only for very few models (see the discussion in Ref. [J|), 
most of unbiased results has been obtained from ad- 
vanced numerical approaches like exact diagonalization 
(ED) [5| , time-dependent density matrix renormalization 
group (tDMRG) [6] or nonequilibrium dynamical mean- 
field theory (nDMFT) ^^l- These approaches allow for 
studying various phenomena, e.g., the nonlinear trans- 
port through interacting nanosystems [8| , the electric- 
field induced breakdown of the Mott insulator [5|, la] , or 
Bloch oscillations [71, llO|. In most cases, theoretical pre- 
dictions for the real-time response can be formulated 
only on the basis of numerical studies. 

It is understood that for a weak driving force, the real- 
time response is determined by the equilibrium correla- 
tion functions as described by the linear-response (LR) 
theory. This theory has recently been extended to ac- 
count for a system, that was previously driven out of 
equilibrium [ll]. Furtheron we focus on the case of the 
charge current I{t) in a one-dimensional (ID) system as 
induced by electric field F{t), whereby the relevant equi- 
librium LR function is the dynamical conductivity cr(a;). 
For a generic (non-integrable) system of interacting (cor- 
related) fermions at finite temperature T > one ex- 
pects within the LR regime the relaxation of current due 
to (Umklapp) scattering processes and hence finite d.c. 
value (Jo = o'(O)- This, in turn, leads to a steady cur- 
rent under a constant driving field F. In contrast, it has 
been recognized that integrable systems in spite of inter- 
action reveal a dissipationless component of the current 
response even at T > as manifested by a finite charge 



stiffness D{T > 0) > [l^- Hence, the qualitative differ- 
ence between both types of systems is expected to remain 
even for the nonlinear transport. 

In this Letter, we address two aspects of the nonlin- 
ear transport in an isolated system of interacting tight- 
binding fermions under a time-dependent driving force: 
(i) How to generalize the LR response approach of generic 
systems at T > to stronger fields and longer steady 
driving? Here we show that the dominating lowest-order 
mechanism beyond the LR regime is the increase of in- 
ternal system energy or Joule-type heating (although we 
are not dealing with relaxation to canonical equilibrium) 
proportional to the square of the electric field. It could 
be accounted for by renormalization of the kinetic en- 
ergy (as the sum rule for a{ui) and, consequently the rel- 
evant (To)- This extension allows one to predict strongly 
nonlinear response without explicit solution of the von 
Neumann or the time-dependent Schrodinger equations. 
(ii) Is there a qualitative difference in nonlinear response 
between integrable and non-integrable systems? Our re- 
sults reveal a clear distinction between both categories 
whereby the response of integrable system to a steady 
field F{t > 0) = const resembles the (damped) Bloch os- 
cillations of noninteracting (NI) fermions, however, with 
nontrivially modified oscillation frequency. 

We investigate an isolated ID system of charged spin- 
less fermions with periodic boundary conditions. The 
system is threaded by a time-dependent magnetic fiux 
(/)(i) and we assume that the flux enters only the kinetic 
energy term Hk of the Hamiltonian H 



H 



Hk + Hi, 



(t) J 



^j+i^j 



h.c. 



(1) 



where = 4>/L is the flux per bond, L is the number of 
sites and Hj is assumed to be flux independent. We put 
furtheron t^ = 1. The time-dependent flux induces an 
electric field F = —0. The charge current operator can 



be written as 



dH 



oil I Y^ f 






h.c 



■} 



(2) 



Before specifying a particular form of the interaction 
term Hj and discussing numerical results it is instruc- 
tive to derive simple equations of motion for the total 
energy E{t) = {H{t)) = Tr[p{t)H{t)], the kinetic energy 
Ek{t) = (Hkit)) and the current I{t) = (J(t)}. Making 
use of the von Neumann equation ip{t) — [H{t),p{t)] one 
can easily find the relations 



E{t) =. ^tTr{[H(t),p{t)]H{t)} + ^{t)Tr 






= 


= LF{t)I{t), 


(3) 


Ek{t) -- 


= i{[H{t),Hkm+LFit)Iit), 


(4) 


lit) -- 


= ^{[H{t),Jm-F{t)^. 


(5) 




FIG. 1: (Color online) Real-time response: a) the kinetic 
energy Ek{t), and b) the renormalized current I{t)/{Ek{t)F) 
obtained for N = 2G,V = 1.4 and W = 1. The electric field 
F is switched on at t = and kept constant for i > 0. 



In the case of NI fermions Hj = 0,H = Hk and commu- 
tators in Eqs. (|4][5]) vanish. Then, for a constant electric 
field F{t) = F equations lead to harmonic oscillations 
with a frequency uib = F, the solution known as Bloch 
oscillations. 

Eq. ([5]) offers also a nontrivial approach to the nonlin- 
ear inverse problem, i.e., for arbitrary interaction term 
Hi one can find the tuning of the electric pulse F{t) 
so that the induced current follows required time profile 
Ja(i). Namely, putting /a(0 ii^ the Ihs. of Eq. (O one 
can solve the equation for F{t). This value of F{t) de- 
termines the magnetic flux that should be used in the 
subsequent time step of the numerical solution of the 
von Neumann equation. Numerical results based on this 
method will be presented elsewhere. 

In the following we study numerically the real-time 
current response within the IT) t — V — W model. 



Hi = V 



'njhj+i+W\njnj+2, 



(6) 



where n,- = 



c]c,, 



whereas V and W are (repulsive) in- 



teractions among particles on the nearest-neighbor and 
next-nearest-neighbor sites, respectively. In the absence 
of the external field the t — V {W — 0) model is integrable 
and within the LR theory 17(0;) has anomalous properties 
[12| . In particular, it can show a dissipationless compo- 
nent even at T > , i.e. a{uj - 0) = D{T)d{uj). W j^ 
breaks the integrability leading to a finite d.c. ao{T > 
0) < 00 [1^. In the following we focus on the metallic 
regime of the half-filled systems with n — N^/L — 1/2 
and V <2th. 

Our aim is to study primarily a generic situation at 
T > to avoid more specific cases emerging from the 
ground state. We perform the numerical evolution of the 
many-fermion wavefunction |^i(i)). The initial |^;(0)) 



should in principle be chosen as eigenstates of i?(0) dis- 
tributed according to the canonical ensemble. Since the 
latter approach is possible only via full ED (typically 
L < 16 for the problem under consideration), we instead 
work using the microcanonical Lanczos method (MCLM) 
[13i] allowing systems up to L = 26. First, we numeri- 
cally generate approximate |\l/i(0)), Z ~ 1, Ng with energy 
E{0) = (iJ(0)) corresponding to the canonical value for 
given T (and L), but as well with a small energy uncer- 
tainty (5^^ = {[H{0)~E{0)]'^). The time evolution |*i(i)) 
is then calculated by step-vise change of (j){t) in small 
time increments 6t <S^ 1 employing at each step Lanc- 
zos basis (typically Nl — 10) generating the evolution 
^;(i — St)) -^ |*;(i)). The described procedure is very 
robust and can be easily tested by changing Ns,6E,St. 
In the following examples we study rather universal (but 
nontrivial) regime of quite elevated T ^ 5 where we use 
iV^ = 10, ^£' = 0.01, St = 0.01. Note that for systems 
considered (L — 26) number of basis states is typically 
Nst ~ 10^ so that SE is still much bigger than the average 
level distance. 

Let us start with the analysis of the non-integrable 
case where we take V = lA and W = 1. For such 
choice of parameters LR cr(w) [13 is rather featureless, 
i.e. broad in lu and, therefore, we may expect that numer- 
ical results are free of artifacts originating from some pe- 
culiar features of cr(aj). Fig.[l]shows the time dependence 
of the kinetic energy Ek{t) (panel a) and the renormalized 
current I{t)/[Ek{t)F] (panel b) for a system driven by a 
constant field F{t > 0) = F. One can see that —Ek(t) 
goes exponentially to zero with the decay rate oc F^ and 
that the ratio I{t)/[Ek{t)F] is almost constant for longer 
t. The dominating nonlinear effect thus consists in the 
increase of E^. Since the latter also represents the sum 
rule, J a{uj)duj oc \Ek\, we formulate and test the follow- 
ing conjecture: I{t) may be well approximated by the LR 



theory extended to account for the increase of Ek (t) , 



Ek{0) 



hnit) 



dt'a{t - t')F{t'), 



(7) 



(8) 



where the conductivity response a{t — t') (evaluated from 
cr(a;)) and consequently /LR(i) are determined within ini- 
tial equilibrium state. Eq. ([7]) alone does not allow one 
to predict the real-time response without solving the 
time-dependent Schrodinger equation, since Ek (t) is still 
needed. Although not fully evident, the accurate scaling 
presented in Fig. [T^ together with Eqs. ([3]) and ([7]) indi- 
cate that the increase of Ek is just proportional to the 
increase of the total energy 



Ek{t)^lE{t)^lLFlER{t). 



(9) 



Therefore, in the long-time regime of a driving with 
F{t) = F one obtains /ER(i) = (JoF Ek{t) / Ek{Q) and, 
consequently, —Ek{t) oc exp(— ai^^i) with a > 0. The 
coefficient 7 is independent of F and can be estimated 
from the initial equilibrium state. 



Ek{t)~Ek{0) 
E{t) - E{0) 



dEkjO) f dE{0) V 
dT \ dT J 



(10) 



The set of equations ([71ITU1) fully determines lER{t)- 
Therefore, similarly to the LR theory, the real-time re- 
sponse can be calculated without explicit solution of the 
time-dependent problem. From the numerical point of 
view, these equations are not more demanding than the 
LR theory, as the only extension consists in the differ- 
ential equation (0). In Fig. [5] we show the accuracy of 
the approximation given by Eq. ([T]). This figure demon- 
strates also the significance of the increase of total energy 
or Joule heating as the dominating nonlinear mechanism. 
Here, the value of 7 has been determined from full diag- 
onalization of a 10-site Hamiltonian. 

We now turn our attention to the real-time response 
of an integrable interacting system, i.e. W = case, 
where the above extension of the LR theory clearly breaks 
down. Now, a relevant reference is the NI y = system, 
where for F(t) = F current and the kinetic energy ex- 
hibit Bloch oscillations with lub — F. It is expected that 
even at T^ > a dissipationless component of current 
will retain the similarity to a NI system. We present in 
the following numerical results for the integrable metal- 
lic case V = 1,W = 0, where Z? (at T — >■ 00) covers 
approximately half of the total spectral weight of (j{uj). 

Inset in Fig. [3] shows /((/)) with <j) = —Ft for the inte- 
grable t-V system driven by a constant electric field F. 
Contrary to the previously discussed W ^ case, one 
observes a strong oscillatory behavior. For large enough 
F > F* ^ 0.5, the frequency of the current oscillations 
oj (evaluated here from the second maximum of I(t)) is 
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FIG. 2: (Color online) I(t) for N = 26, V = 1.4 and W = 1. 
Arrows indicate the instants of time, when the electric field is 
switched on (all panels) and off (upper left panel). Full lines 
represent results obtained via unrestricted numerical solution, 
whereas dotted and dashed lines show approximate results of 
the LR theory, 7LR(i), and the extended LR approach, /ER(i)i 
respectively. 
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FIG. 3: (Color online) Current-oscillation frequency lj vs. 
field F normalized to ujb of NI fermions, evaluated for A^ — 26, 
V ^ 1,W = 0. Inset: / vs. fiux \4>\ for different F/tt = 
0.02,0.04,0.08,0.16,0.4. 
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FIG. 4: (Color online) Relation between current I{t) and ki- 
netic energy density Ek(t)/L for parameters as in Fig. [3l for 
different fields F/tt = 0.02 — 0.4. Shown is also the results for 
NI fermions. The point and the arrow mark the initial state 
and the direction of evolution, respectively. 



approximately the same as ujb of NI fermions. Simi- 
lar prediction has been obtained within nDMFT (see the 
discussion of a metallic phase in Ref. [7|). However, in 
the regime of a low electric fields F < F* , the corre- 
lated system oscillates faster than ujb ■ Within the inves- 
tigated t-V model, ui/ujb is found to increase logarithmi- 
cally when the electric field decreases (see Fig. [3]) . The 
effect is present also for other V < 2th although evidently 
the logarithmic regime disappears for y — > where the 
NI universality sets in, i.e. in this case the crossover also 
moves to F* ^ 0. 

At present we are unable to provide a physical explana- 
tion for the enhancement of uj/lob or to predict whether 
the logarithmic scaling breaks down for very low F. At 
least it is evident that F > F* response is similar to NI 
fermions in that both I{t) and Ek{t) — E{t) change sign 
during the evolution (see Fig. HJ. On the contrary, the 
kinetic energy remains negative in the weak-field regime. 

Still, for sufhciently short times of driving i > the 
system's state is almost the initial one, i.e. p{t) ~ /^(O), 
and the changes of J(i) and Ek{t) come from their explicit 
flux dependencies in Eqs. ([TJS]). Then, one can easily 
find from Eqs. dUE]) that [Ek{t)Y + [LI[t)Y ~ [EkiG)]^ , 
independently of the form of interaction term Hj. It 
explains, why the initial response as well as the strong- 
field response of a correlated system are the same as the 



response of NI fermions [ij]. Moreover, assuming that 



the Ek{t) and I{t) do not cross the outer circle in Fig.|4] 
(what holds true in all investigated cases), one finds that 
the magnitude of current is bounded from above by the 
initial kinetic energy. 

In conclusion, our study reveals several novel features 
of isolated interacting fermionic systems under constant 
or time-dependent driving force. For a generic system 
of tight -binding electrons we present a simple extension 
of the LR theory taking into account the change of ki- 
netic energy. This approach gives a very satisfactory de- 
scription of numerical results, at least for high enough 
T . Here, the basic condition is that in spite of driving, 
the system at all times satisfies the quasi-equilibrium re- 
lation between the kinetic energy Ek{t) and the total 
energy E{t). It seems plausible that the necessary condi- 
tion for such a development is (fast enough) relaxation of 
current, here due to the Umklapp processes. It remains 
to be investigated whether a similar behavior applies to 
interacting systems beyond the tight-binding description 
(when e.g. the sum rule for a{uj) is not directly Ek) or 
in low T regime. In the latter case the linear dependence 
Ek{t) ex E{t) should probably be replaced with a more 
general one Ek{t) — Ek[E{t)]. 

On the other hand, integrable systems show strikingly 
different development. To first approximation the cur- 
rent I{t) under constant driving reveals oscillations with 
the dominant dependence only on the fiux ^(i). While 



at large F > F* this has clear connection to the Bloch 
oscillations of tight-binding electrons with characteristic 
Lu ^ UJB, the effective uj is increasing for lower F < F* 
while at the same time oscillations are becoming damped 
(in fact the damping is found to be strongest at F ~ F*). 
One part of this phenomena, in particular the damp- 
ing, can be attributed to the breaking of integr ability 
by introducing time-dependent 4>{t). Oscillations them- 
selves in the low-field regime are the signature of a co- 
herent behavior which has to be intimately related to 
finite charge stiffness D{T > 0) > since they disap- 
pear outside the metallic regime, e.g. for V > 2th at 
half-filing n = 1/2. Alternatively, if the field is switched 
off in this regime before reaching Ek ^ then the cur- 
rent relaxes to a finite I{t = cxd) > being a direct 
consequence of D(T > 0) > 0. In any case, the ob- 
served phenomena reveal that even in weak fields the 
response of (near)integrable systems is far from the one 
expected from the equilibrium LR theory, but at the same 
time quite universal representing an open challenge for a 
proper explanation and also possible experimental obser- 
vation. 
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